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ABSTRACT 



More than half of stars reside in binary or multiple star systems and many 
planets have been found in binary systems. From theoretical point of view, 
^ ■ however, whether or not the planetary formation proceeds in a binary system is 

^ ■ a very complex problem, because secular perturbation from the companion star 

<3\ • can easily stir up the eccentricity of the planetesimals and cause high-velocity, 

^ ■ destructive collisions between planetesimals. Early stage of planetary formation 

■ process in binary systems has been studied by restricted three-body approach 

^ ■ with gas drag and it is commonly accepted that accretion of planetesimals can 

proceed due to orbital phasing by gas drag. However, the gas drag becomes 



less effective as the planetesimals become massive. Therefore it is still uncertain 
whether the collision velocity remains small and planetary accretion can proceed, 
once the planetesimals become massive. We performed A''-body simulations of 
planetary formation in binary systems starting from massive planetesimals whose 
size is about 100-500 km. We found that the eccentricity vectors of planetesimals 
quickly converge to the forced eccentricity due to the coupling of the perturbation 
of the companion and the mutual interaction of planetesimals if the initial disk 
model is sufficiently wide in radial distribution. This convergence decreases the 
collision velocity and as a result accretion can proceed much in the same way 
as in isolated systems. The basic processes of the planetary formation, such as 
runaway growth and oligarchic growth and final configuration of the protoplanets 
are essentially the same in binary systems and single star systems, at least in the 
late stage where the effect of gas drag is small. 
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Subject headings: binaries: close — planetary systems: formation — methods: 
n-body simulations 

1. Introduction 

As of March 2007, 215 candidates of extra-solar planets have been found and at least 
30 of them are in binary or multiple star systems (Raghavan et al. 2006). Table [1] shows the 
candidates of the close binary or multiple star system which has the planets. In this table, 7 
Cephei is an example of close binary systems in which we are interested. According to Hatzes 
et al. (2003), the semi-major axis and eccentricity of the companion star of 7 Cephei are 
18.5 AU, and 0.36, respectively. Planetary formation process in such a close binary system is 
the main target of this paper. The frequency of planets in binary systems is not significantly 
different from that for single star systems (Desidera and Barbieri 2007). It is very important 
to investigate the formation process in binary system because more than half of stars reside 
in binary or multiple star systems. 

Many authors have investigated planetary accretion process using A'"-body simulation 
(e.g., Aarseth, Lin, and Palmer 1993; Kokubo and Ida 1996, 1998, 2000, 2002). In all of 
these simulations, the evolution of protoplanetary disks around isolated stars was studied. 
On the other hand, Quintana et al. (2002, 2006) have investigated the planet formation in 
binary system from protoplanets by A''-body simulations. They studied the late stage of the 
planet formation process (from protoplanets to planets). 

Marzari and Scholl (2000) and Thebault et al. (2004, 2006) investigated the distribution 
of the collision velocity between planetesimals with and without gas drag in binary systems 
by restricted three-body approach. They found that orbital phasing occurred and collision 
velocity remained small due to the coupling of gas drag effect and secular perturbation, 
and they concluded that planetary accretion could proceed. As the planetesimals become 
massive, however, the gas drag becomes ineffective and mutual gravitational effect between 
planetesimals becomes important. In such a condition, whether the collision velocity remains 
small is not clear. As stated above, the final stage, from protoplanets to planets, has been 
studied by A^-body simulations, but there has been no A^-body work on the intermediate 
stage from massive planetesimals whose size is 100 km - 500 km to protoplanets. In this 
paper, we focus on this intermediate stage. 

In a binary system, the orbits of planetesimals change due to perturbation from the 
companion. The most important term of the perturbation in our simulation region is secular 
perturbation, if the orbit of the companion is eccentric. The effect of secular perturbation 
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is that eccentricity vector moves on the "perturbation circle". Thus, if interaction between 
planetesimals and gas drag are neglected, planetesimals with different semi-major axis move 
on the perturbation circles on different frequencies and phases, and therefore gain high 
relative velocity. This is the reason why the collision velocity becomes high in previous 
works without gas drag or self gravity (Marzari and Scholl 2000; Thebault et al. 2004, 
2006). If the collision velocity becomes high when the planetesimals become massive and 
thus, the effect of gas drag becomes small, accretion process might halt since collisions might 
be destructive. 

It is a very important question whether destructive collisions occur when the gravita- 
tional interaction between planetesimals is taken into account because it determines whether 
the accretion process can continue after gas drag becomes ineffective. Ito and Tanikawa 
(2001) studied the evolution of the orbital elements of protoplanets under the perturbation 
of Jupiter and found that the eccentricities of the protoplanets aligned with each other, 
resulting in the evolution was very similar to that without the influence of Jupiter. This 
alignment is due to gravitational interaction between the protoplanets. If this kind of align- 
ment also occur for planetesimals in binary system, the collision velocity can become smaller 
than the value predicted by Marzari and Scholl (2000). 

In this paper, we report the results of A^-body simulations of protoplanet formation from 
massive (3 x 10^"^ — 1.4 x 10^^ g) planetesimals. We start simulations with the planetesimal 
disk in which all planetesimals have the same mass. The gravitational interaction between 
planetesimals is included and gas drag is neglected. We summarize theory of secular per- 
turbation and discuss what initial distribution should be used in section 2. The numerical 
scheme and initial conditions are described in section 3. We also discuss the vahdity of the 
"same mass" assumption of planetesimals in section 3. The results of A^-body simulations 
are presented in section 4. In section 5, we sum up. 

2. Theoretical Prepciration 

2.1. Brief description of secular perturbation 

If the mutual interaction of planetesimals and gas drag are negligible, the orbital evolu- 
tion of a planetesimal in a binary system can be described by restricted three-body approxi- 
mation, and the time cvohition of eccentricity and longitude of pcriccntcr of the planetesimal 
arc given by secular perturbation theory. The secular evolution of the h and k variables of 
the planetesimal are expressed as (Heppenheimer 1978; Whitmire et al. 1998; Thebault et 
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Table 1: Known candidates of close systems 







a sm i 
(AU) 


Projected 


M sin i 
(Mj) 






Name 


Component 


Separation 
(AU) 


eccentricity 


References 


HD 041004 


B 




22 






1, 2, 3, 4 




C 


0.016 




18.4 


0.08 


3,4 




b 


1.31 




2.3 


0.39 




7 Cephei 


B 


20.3 






0.39 


1, 5, 6, 7, 8 




b 


2.03 




1.59 


0.2 




GJ 893 


B 
C 
b 


0.3 


2248 
18 


2.9 




9, 10 



Note. — These parameters are from Raghavan et al. (2006) and its references. Column component lists 
companions (B, C, D,...) and planets (b, c, d, ...). 

References. — (l)Eggenberger et al. 2004; (2)Sce 1896; (3) Zucker et al. 2003; (4) Zucker et al. 2004; (5) 
Mason et al. 2001; (6) Campbell et al. 1988; (7) Griffin et al. 2002; (8) Hatzes et al. 2003; (9) Zacharias et 
al. 2004; (10) Wilson 1953 

al. 2006) 

h{t) — epSm{At -\- wq) , (1) 

k{t) — epCOs{At -\- voq) -\- tf , (2) 

where 

h{t) — esin(ro) , (3) 

k{t) — ecos(ro) , (4) 

and e and w are the eccentricity and the longitude of pcriccnter of the planetesimal. Here, 
^, Cp, -Too, and ej are constants which depend on the strength of the secular perturbation. 
We take the k axis as direction of the eccentricity vector of the companion. In other words, 
the eccentricity vector of the companion is (e^, 0). The semi-major axis of the companion is 
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qb and its mass is Mb- According to the secular perturbation theory, e/ is given by 

_ 5 a cb , . 

This 6/ is forced eccentricity induced by the companion star. Here, a is the semi- major axis 
of the planetesimal. The angular velocity of rotation on k-h plane. A, is given by 

3 1 a^/^ 



2 (l-e|)3/2 



a 



We use system of units in which the solar mass, 1 AU and 1 year are all unity. Remaining two 
constants, zuq and Cp, are determined from the initial orbital elements of the planetesimal. 
We call the vectors e = {k,h), e/ = (e/, 0) the eccentricity vector and the forced eccentricity 
vector, respectively. On k-h plane, the eccentricity vector of a planetesimal moves on a circle 
centered at e/. 



2.2. Initial eccentricity of planetesimals 

We consider two types of initial conditions for planetesimals. In the first one, the initial 
distribution of eccentricity vectors of planetesimals is centered at the forced eccentricity 
vector. The distributions of e — e/ and inclination i are both given by Rayleigh distribution 
with dispersions (|e — e/p)^/^ = 2(i^)^/^ = 0.02. We call this model "forced" model. In the 

second model, the distribution of eccentricity vector is centered at the origin of k-h plane. 
The distributions of eccentricity e and inclination i are also given by Rayleigh distribution 
with dispersions (e^)!/^ = 2(^2)1/2 = o.02. This is the same distribution as those used in 
previous works for iV-body simulation of planetary formation in isolated systems (Kokubo 
and Ida 2002). We call this model "circular" model. 

We argue that the "forced" model is more suitable for simulation in binary system than 
the circular model for the following reasons. Consider a narrow region of a planetesimal disk 
such as the region of 0.95 AU < a < 1.05 AU. If the eccentricity vectors of planetesimals in 
this region do not distribute around the forced eccentricity vector, they will rotate around the 
position of the forced eccentricity vector. The planetesimals in a narrow region would rotate 
together because the angular velocity A is similar for planetesimals with similar values of the 
semi-major axis a. In addition, they align due to secular interactions between planetesimals 
(sec, e.g., Ito and Tanikawa 2001), Because of this collective motion, the relative velocity 
between planetesimals is kept small in this narrow " ring" . 

However, if we consider a wider region, it would behave as collection of many narrow 
rings. If we ignore interactions between the rings, they would rotate on their own angular 
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velocities A. Since neighboring rings have shghtly different values of A, they would soon 
physically collide with each other, resulting in damping of "free" eccentricity. Gravitational 
interaction between rings would also dump relative difference of eccentricity vectors. In other 
words, planetesimals in the circular model would first relax to the forced model. We argue 
that the equilibrium state is more suitable for the initial conditions. This is why we consider 
"forced" model. 

As we mentioned in introduction, we study the late stage of the formation process of 
protoplanets. The earlier phase have been studied by Marzari and Scholl (2000) and others. 
Thus, it might seem reasonable to set the initial eccentricity and longitude of pericentcr of 
planetesimals to the equilibrium value of the orbital phasing in Marzari and Scholl (2000) 
instead of the circular model. However, We found it is a bit difficult to use their results, 
because they adopt circular gas disk model. This circular gas disk in a binary system with 
eccentric companion seems a bit unnatural. Due to the perturbation from the companion 
star, the gas disk may be twisted. We could not find previous works which directly studied 
the equilibrium state of gas disk under the secular perturbation of companion, but recent 
work by Papaloizou (2005) seems to imply that eccentric disk can be long-lived even if there 
is no companion. So it seems likely that gas disk is not circular when eccentric companion 
exists. This may change gas drag effect to the planetesimals and might affects the direction 
of the orbital phasing. This is the reason why we do not use the Marzari's results. Even if 
the orbital elements of them is correct, our circular model is very close to their results and 
it would soon relax to the "forced" model. So we expect that it makes no significant change 
to our results. Clearly, more studies on the dynamics of gas disk is necessary. 



2.3. Collision velocity 

As noted in section 2.1, we adopt the Rayleigh distribution with dispersions (e^)^/^ = 
2(^2)1/2 ^ 0.02 for the "circular" model (Ida and Makino 1992), and (|e-e/H V2 ^ 2(^2) V2 ^ 
0.02 for the "forced" model (see Table 2 ). These conditions imply collision velocity is 
initially Vcoi — 500 — 1000 m s~^ at 1 AU for the initial distribution. When the effect of 
gas drag is taken into account, the equilibrium eccentricity and inclination are estimated as 
(e2)i/2 = 2(^2)1/2 = 0.0042 (Kokubo and Ida 2000). This value is probably more reasonable. 
In this paper, however, we use the value 0.02 because we neglect gas drag force in our 
calculation and the value of the eccentricity and inclination will soon relax to the value of 
gas-free case. Our focus is not on the quantitative analysis of the collision velocity but on 
the qualitative analysis of the coupling of secular perturbation and gravitational interactions 
between planetesimals. To determine the realistic value of the coUision velocity, we should 
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include gas drag force in our calculation. 

3. Method of Calculation 

3.1. Initial conditions 

The model parameters for the companion star are summarized in table O Mass of the 
primary star is 1 Mq for all models. We adopt three models for companion. Model alpha 
corresponds to a Cen system. The mass ratio is chosen to be the same as that of a Cen 
system, (1.1:0.9). For the other two models, we use a somewhat smaller semi-major axis 
than those of observed binary systems in which the planets are found. We choose this value 
to study the case in which secular perturbation is strong. For circular models, we use an 
axisymmetric surface mass density distribution for planetesimal disk whose surface mass 
density is given by 

^soud = ^oij^)-'^' g cm-2 , (7) 

where a is distance from the primary star, and Sq is the reference surface density at 1 AU 
(We adopt Sq = 10 g cm~^). Table [3] shows model parameters. Here, is the number of 
planetesimals. For radial distribution, we use two disk models, "wide disk" models (models 
0-5) and "narrow disk" models (models 6-8). In the "wide disk" models, we set inner 
and outer cutoff radii to be 0.5 AU and 1.5 AU, respectively. In the "narrow disk" models, 
inner and outer cutoff radii are 0.95 AU and 1.05 AU, respectively. Planetesimals have equal 
mass in all models. The number of planetesimals is 10000 in the wide disk models. In two 
narrow disk models (model 6, 7), Number of planetesimals is 5000. In one of narrow disk 
models (model 8), the number of planetesimals is 975. This is a "cutoff" model which has 
the planetesimals of the same mass as in the wide disk models. The density of planetesimals 
is 2 g cm~^. We increase their radii by a factor 5 to accelerate accretion process (Kokubo 
and Ida 1996). At first, all planetesimals have same mass (~ 1.44 x 10^^ g in models 0-5 
and 8 and 2.88 x 10^^ g in model 6 and 7). 

One critical question is if our " same mass" setup can be really regarded as description of 
the later stage of planetary formation. Recent theoretical and numerical works of planetary 
accretion (for example Inaba et al. 2001, Rafikov 2003) seem to suggest that the "orderly" 
phase of growth does not exist, and runaway growth starts from much smaller mass than 
that of our setup. Our setup, therefore, is not realistic. However, statistical calculation 
by Inaba et al. (2001) has shown that once massive planetesimals (M ~ 10^^ g) formed, 
planetesimals with mass less than 10^^ g become dynamically unimportant. So even though 
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our initial condition is oversimplified, it might still give qualitatively valid description of the 
later phase. 

Two planetesimals are considered to coUide when their distance becomes less than the 
sum of their radii. We assume the perfect accretion under which planetesimals always accrete 
when they coUide. Whether or not this assumption is good depends on the distribution of 
the collision velocities. The escape velocity of our planetesimals is 200 — 600 m s~^ and is 

of the same order as the initial collision velocity. Furthermore, if gas effect is included, the 
eccentricity dispersion might be less than our estimate and the collision velocity becomes 
smaller. So we believe the assumption of perfect accretion is valid. 

As stated above, we increase the radii of planetesimals ^fold (here, / = 5) to save 
calculation time. This acceleration of accretion may affect time evolution of planetesimals, 
especially the time evolution of orbital elements on k-h plane because the increase of radii 
changes timescale of accretion but does not affect secular perturbation of companion star. 



Table 2: Model parameters for the companion star 



semi-maior axis . . 

Name Mass (AU) eccentricity 

nocomp 

e25 IMq 16 0.25 

e50 IMo 16 0.5 

alpha O.82M0 23.4 0.52 



3.2. Integration scheme 

We use the fourth-order Hermite scheme (Makino and Aarseth 1992) with hierarchical 
timesteps (Makino 1991) improved for planetary systems (Kokubo et al. 1998 ) for numerical 
integration of planetesimals and companion star. The equation of motion for planetesimals 
is given by 

a, = - 5^ Gm,^ - G{m, + M,)^ - GmA + ^) , (8) 

i^j ij ip ic ' pc 

where Mp, Mc, rip, and r^c are mass of the primary star, mass of companion star, position of 
planetesimal relative to the primary star and that relative to the companion star. We use 
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Table 3: Initial conditions 



Model 




Width of disk 
(AU) 


companion 


initial orbit 





10000 


0.5-1.5 


nocomp 




1 


10000 


0.5-1.5 


e25 


forced 


2 


10000 


0.5-1.5 


e50 


forced 


3 


10000 


0.5-1.5 


e25 


circular 


4 


10000 


0.5-1.5 


e50 


circular 


5 


10000 


0.5-1.5 


alpha 


circular 


6 


5000 


0.95-1.05 


e25 


circular 


7 


5000 


0.95-1.05 


alpha 


circular 


8 


975 


0.95-1.05 


e25 


circular 



position of the primary star as origin of the coordinate for planetesimals. The motion of the 
primary star due to the gravitational forces of planetesimals is neglected. 

Most expensive part of the numerical integration is calculation of mutual gravitational 
interaction between planetesimals. We use GRAPE-6 (Makino et al. 2003) and GRAPE-6A 
(Fukushige et al. 2005) to calculate the gravitational interaction between planetesimals. We 
also integrate orbit of the companion star using the fourth-order Hermite scheme. For both 
the planetesimals and the companion stars, we use the standard timestep criterion (Aarseth 
1985) 



At = Jr] 



\a 


.||a(2) 


+ 


a 


2 


\a\ 


|a(3)| 


+ 


|a(2)|2 



(9) 



Since the orbital period of the companion is long, the accuracy of its orbit is more than 
enough. 

4. Results 

4.1. Planetary accretion in binary systems 

Figure [T] shows the time evolution of planetesimals of model 1 (left) and model (right) 
on the a-e plane. We integrate the system for 5 x 10^ years. We can see that protoplanets 
grow in a very similar way in these two models. This behavior is essentially the same for all 
models. 
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Figure [2] shows cumulative mass distribution of planetesimals of the region 0.9AU < a < 
1.1 AU. The time evolutions of all of these models are very similar. In all cases, the mass 
distributions first relax to the power-law distribution with power index dlognd dXogm ~ 
— 1.5 (top panels) where Uc is the cumulative number of planetesimals. This power-law 
distribution is characteristic of the runaway growth (Makino et al. 1998, Kokubo and Ida 
2000). 

As time goes on, the mass distributions become bimodal: planetesimals with mass 
M ~ 1 — 10 X 10^^ g and large protoplanets with mass M ~ 10^"^ g. This bimodal mass 
distribution is the result of the runaway and the oligarchic growth of protoplanets (Kokubo 
and Ida 1996,1998). The time evolution of the mass distribution of binary systems is almost 
the same as that in single star systems. It means that secular perturbation does not change 
the process of the protoplanet formation. Runaway growth and oligarchic growth also occur 
in binary systems in a similar way as in the single star system. 

In figure [3] - [6l we show the evolution of the mass of most massive planetesimals ([3] and 
H]) and average mass of planetesimals ([S] and E]) • We can see that the evolution is rather 
similar and there is no systematic tendency due to the presence of the companion. 

The evolution of the surface mass density of model 1 (e = 0.25, forced) and that of model 
(no companion) are shown in Figure [71 We found that density profile of planetesimals in 
a binary system is virtually the same as that in a single star system. 

From secular perturbation theory, time evolution of semi-major axis of a planetesimal 
is given by 

da 2 dR , , 

dt na oe 

where e and R are mean longitude at epoch and disturbing function, respectively. In practice, 
variation of e can usually be neglected since it is a small effect. Thus, mass migration does 
not occur since secular perturbation hardly change the semi-major axis of planetesimals. 
Our results are consistent with this theoretical expectation. 



4.2. Time evolution of eccentricity 

The time evolutions of the distribution of the planetesimals on k-h plane of model 1 
(forced model of e = 0.25), 3 (circular model of e = 0.25), 8 (cutoff model of e = 0.25), 4, 
and 5 (circular models) are shown in Figure [H 

The top three panels show models 1, 3, 8, all with an e25 companion. In model 1 (top 
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panel), planetesimals are initially distributed around the forced eccentricity, and this does not 
change in time. In model 3, planetesimals are initially distributed around zero eccentricity, 
but this initial distribution is replaced by a distribution centered at the forced value. The 
behavior of planetesimals in model 8 is very different. The distribution after 30,000 years 
is not centered at the forced value and keeps rotating around the forced eccentricity vector 
(see figure [9] and description below). This difference between narrow disk models and wide 
disk models will be discussed in more details in section 4.2.1. The bottom two panels of 
figure M show that the eccentricity vectors of the planetesimals of other circular models with 
wide distribution also converge to the forced eccentricity vector. We will discuss the time 
evolution of proper eccentricity in section 4.2.2. 



4.2.1. Difference between wide and narrow models 

We performed simulations of the narrow disk models (models 6, 7, 8) to see the differ- 
ence between the wide and narrow distributions. We use two types of narrow disk models. 
One has the planetesimals of the same mass as in wide disk models, and thus the number 
of planetesimals is small (A^ = 975, model 8). The other has a large number of smaller 
planetesimals (see Table 2). 

Figure M shows the time evolution of eccentricity. In wide disk models, the oscillation 
of eccentricity is damped quickly due to the convergence to the forced eccentricity vector 
as shown in figure [81 In narrow disk models, however, the eccentricity liberates around the 
forced eccentricity and does not converge to the forced eccentricity. 

This difference can be understood as follows. As we discussed in section 2, the planetes- 
imals rotate around the forced eccentricity vector on k-h plane and the angular velocity A is 
a function of semi-major axis of planetesimals. In the case of narrow disk models, the range 
of A is small, and it is possible that all planetesimals synchronize due to mutual gravitational 
interaction. Thus in narrow disk models planetesimals move collectively on k-h plane. In the 
case of the wide disk models, however, such collective motion is not allowed since the range 
of A is too large. If planetesimals with different values of a (therefore, ^4) rotate around 
the its own values of forced eccentricity on their own timescales, collision is enhanced and 
eccentricity will be damped to forced values. Thus, oscillation of eccentricity damps quickly 
in wide models. 

We performed a simulation of narrow disk models to see how this behavior is affected 
by the mass of planetesimals. The bottom panel of Figure [9] shows the time evolution of 
eccentricity of model 8. The behavior of model 8 is almost the same as that of model 5. 
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Thus, the number of planetesimals or mass of each planetesimal do not change the resuh. 
In narrow disk models, the convergence to the forced eccentricity vector never occur. 

4-2.2. Time evolution of proper eccentricity 

Time evolution of proper eccentricity, {i.e., distance from the forced eccentricity in k- 
h plane) is shown in figure [TOl The proper eccentricity is averaged for planetesimal with 
semi-major axis between 0.95 AU < a < 1.05 AU. In circular models (model 3, 4, and 
5), the proper eccentricity quickly decreases through coUisional damping and gravitational 
relaxation. This means that the eccentricity vector of the planetesimals converge to the 
forced eccentricity vector. In forced models (model 1, 2) and single star model (model 0), 
on the other hand, the proper eccentricity does not change significantly in early stage of the 
simulation (before 40,000 years). The convergence in model 5 {a Cen model) is slower than 
that in the other two circular models and before convergence the planetary accretion has 
almost finished. This is because of the small angular velocity ^ of a Cen system. However, 
if we adopt real radii, the convergence would probably occur before the accretion completes. 
The time evolution of proper eccentricity after the convergence is similar to that in the forced 
models and the single star model. The proper eccentricity increases in the late stage due to 
the viscous stirring. 

Figure [11] shows time evolution of the mass weighted average of proper eccentricity for 
the same radial range as in figure [TOl In the early stage of the simulation, its behavior is 
similar to that of the simple average in figure [TOl However, in the late stage, the increase 
is slow because of dynamical friction to seeds of protoplanet. In single star systems, the 
eccentricity of seeds of protoplanet becomes small {i.e., the eccentricity vectors converge 
to the origin of k-h plane). On the other hand, in binary systems, the eccentricity vectors 
converge to the forced eccentricity and the evolution of proper eccentricity is similar to the 
evolution of eccentricity in the single star system. This means that the orbit of planetesimals 
does not become circular as in single star system but become eccentric in binary systems 
and their pericenter aligns to the pericenter of companion star. 

The distance between planetesimals on k-h plane determines the collision velocity. Thus, 
the convergence reduces the collision velocity between planetesimals. In model alpha, for 
example, the maximum distance between planetesimals could becomes 0.08 by secular oscil- 
lation if the mutual gravitational effect is neglected. It corresponds to the collision velocity 
of 3000-4000 m s~^ at 1 AU. Collision with this velocity is destructive even for the mass of 
the planetesimal of about 1.4 x 10^^ g (its escape velocity is about 600 m s~^). On the other 
hand. The collision velocity is reduced to about 500-1000 m s~^ with this convergence. As 
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we mentioned above, gas drag affects tfie eccentricity dispersion of the planetesimals {i.e., 
gas drag can not be negligible in this sense). The simulations in which gas drag and mutual 
gravitational effects are taken into account is required to determine the precise value of the 
collision velocity. 

5. Conclusion and Discussion 

We have performed the simulations of formation of protoplanets from massive plan- 
etesimals (with radii of 100km- 500km) in close binary systems, for which gas drag effect 
is negligible and the coupling of the secular perturbation and gravitational interaction is 
important. 

We found that eccentricity vectors of planetesimals quickly converge to the forced ec- 
centricity if the initial disk model is sufficiently wide in radius and secular perturbation is 
sufficiently strong. The eccentricity vectors of planetesimals which have heavier mass than 
the other planetesimals converge to the forced eccentricity more strongly. This convergence 
results in orbital phasing and it reduces value of collision velocity to less than the value 
predicted in studies with restricted three-body approach. Runaway growth and oligarchic 
growth also occur in binary system much in the same way as isolated star system at least 
in the late stage of the formation. The final configuration of protoplanets is not different 
between close binary systems and isolated star systems. 

Our simulations, however, have following limitations. (1) We underestimated effect 
of secular perturbation. The /-fold change of radius increases frequency of collision and 
accelerates planetary accretion. This acceleration of the formation process causes relative 
under-estimate of secular perturbation. We plan to perform the simulation with real radii 
of planetesimals to sec if this effect would make any difference. (2) We investigated only 
three binary systems. More simulations of various binary systems are required to determine 
quantitative relationship between the convergence of eccentricity vectors and strength of 
secular perturbation. (3) We neglected effect of gas drag. This is okay for the late stage 
which we studied. In earlier stage, however, gas drag is clearly important. We plan to 
calculate the equilibrium state of gas disk under the secular perturbations and perform the 
simulations which include the gas effect in future works. 
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Fig. 1. — Distribution of planetesimals on the a-e plane for model 1 (left) and model 
(right). The circles represent planetesimals. The filled circle represent protoplanets with 
mass more than 100 times the initial mass. 
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Fig. 2. — Cumulative mass distribution of planetesimals. Left panels show the results of 
the forced models (model 1, 2) and right panels show the results of the circular models 
(model 3, 4). In all panels, the results of isolated model (model 0) are also shown. Times 
are 5 x 10'', 1 x 10^, 3 x 10^, 5 x 10^ years from top to bottom. The dotted hues at the top 
panels indicate the slope of -1.5. 



-18- 



1000 
900 
800 
700 
600 
500 
400 
300 
200 
100 




model 
model 1 
model 3 




I 

J 



50000 100000 150000 200000 250000 300000 

time (yr) 



Fig. 3. — The time evolution of maximum mass of planetesimals between 0.9AU < a < 
I.IAU. Solid, dashed, and dotted curves show the result of the model 0, 1, and 3, respectively. 
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Fig. 4. — Same as figure [3] but for models 0, 2, and 4. 
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Fig. 5. — Same as figure [3] but the average mass is shown. 
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Fig. 6. — Same as figure [5] but for models 0, 2, 4. 
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Fig. 7. — The surface mass density of planetesimals for model 1 (solid) and model (dotted) 
at t = 1 X 10^, 1 X 10^, 2 X 10^, 4 x 10^ years, from top to bottom. 
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Fig. 8. — The distribution of planetesimals of model 1, 3, 8, 4, 5 (top to bottom) on k-h 
plane. Left panels show initial distributions and right panels show evolved state (3 x 10^ 
years, for models 1, 3, and 8 and 2 x 10^ and 8 x 10^ years for models 4 and 5) We plotted 
planetesimals between 0.95 AU < a < 1.05 AU. The center of circles is the forced eccentricity 
for a = 1 AU and radii of circles are 0.01, 0.03, 0.05 respectively. The open and filled circles 
represent the planetesimals with mass 10 times and 100 times the initial value. 
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Fig. 9. — The time evolution of the average eccentricity in narrow disk models (models 6 - 
8) compared with those in wide disk models (models 3, 5, and 3 from top to bottom). For 
wide disk models, the eccentricity is averaged for planetesimals in the range of 0.95 AU < 
a < 1.05 AU. The sohd curves show results of wide disk models and dashed curves show 
results of narrow disk models. 
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Fig. 10. — The time evolution of average of proper eccentricity, the proper eccentricity is 
averaged for planetesimals with 0.95 AU < a < 1.05 AU. Sohd, long-dashed, short-dashed, 
dotted, and upper dot-dashed curves show the results of models 1-5, respectively. Lower 
dot-dashed curve shows the result of model 0. 
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Fig. 11. — The same as figure ITOl but mass- weighted averages are shown. 



